/**************************************************************************//**
 * @file    Helpers.c
 * @author  İsmail SEZEN (email: sezenismail@gmail.com)
 * @date    02-08-2009
 *
 * @note This is part of siklon project. (Siklon.c)
 *
 * <b>Description</b>
 *
 * General helper routines for siklon project.\n
 *
 * Copyright (C) 2009 İsmail SEZEN
 *
 * <b>License</b>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 *****************************************************************************/

/* Include Directives *///{
#include <stdio.h>
#include <stdlib.h>
#include <logger.h>
#include "Helpers.h"
#include "siklon_internal.h"
#include "array_helper.h"
#include "Parabola.h"
//}

/*Macros*///{
/*! @brief Resets the axis values to default initial values because Min_Search_Parabolic_Interpolation function changes them. */
#define RESET_AXIS(a) a[0] = -2.5; a[1] = 0.0; a[2] = 2.5;
/*! @brief Checks an array elements satisfies a condition like x0 < x1 and x1 < x2. */
#define COND1(a) ((a[0]<a[1]) && (a[1]<a[2]))
/*! @brief Checks first 3 elements of an array are equal or not. */
#define ARE_EQUAL(a) ((a[0]==a[1]) && (a[0] == a[2]))
/*! @brief Takes averages of an array's first and third element. */
#define AVE(a) ((a[0]+a[2])/2)
/*! @brief Gets the value of an 2D array's i. and j. element by using it's pointer */
#define VAL(h,i,j) (*(h + (j*144+i)))
/*! @brief this is a helper macro for for loop. */
#define loop(i,x,y) for((i)=(x); (i)<=(y); (i)++ )
/*! @brief Returns the absolute value of i. */
#define ABS(i) (i)<0 ? -(i) : (i)
//}

/*! @brief This is used for minimizing process*/
double tolerance = 1.0e-10;

/*! @name Extern Functions defined in another file
 *  These functions are defined in other files.
 *  @sa min_search_parabolic_interpolation.c
 *  @sa biparabolic_interpolation.c */
//@{//{
extern int Min_Search_Parabolic_Interpolation(double (*f)(double), double* x0, double* fx0, double tolerance );
extern double Find_Zxy_Biparabolic_Interpolation(double x, double y, double *xi, double *yi, double *zij);
//}//@}

/*! @name Functions related to parabolic calculations and Min/Max detection */
//@{//{

/*! @brief Creates a 3D data array to hold the geopotantial height values.
 *  @param[in] level_count Level count
 *  @param[in] lat_count Latitude Count
 *  @param[in] lon_count Longtitude Count
 *  @return 3D Data Array */
short*** createDataArray(int level_count, int lat_count, int lon_count) {
	return arr_3d_short(level_count, lat_count, lon_count);
}

/*! @brief Tests given value is minimum or not.
 *
 *  Tests the geopotantial height at the given point is less than or equal to
 *  each of the 8 neighboring height-values.
 *  @param[in] *h Pointer to the geoptantial height value in the 3D array
 *  @return Returns an iteger value shows SUCCESS or FAIL.
 *  @retval  0 Height value is NOT minimum (FAIL)
 *  @retval -1 Height value is minimum (SUCCESS) */
int hgt_is_min(short *h){
	register int i,j;
	for(i = -1;i <= 1;i++)
		for(j = -1; j <= 1;j++)
			if( VAL(h,i,j) < *h ) return 0;//Middle Value is not minimum
	return -1;//Yeah! That's it!
}

/*! @brief Search the minimum value and location
 *
 *  This function searches for a minimum point value and location on a single axis.(e.g. x or y)
 *  Internally, uses the \ref Min_Search_Parabolic_Interpolation function to calculate the minimum point location and value.
 *  \ref Min_Search_Parabolic_Interpolation function gives the result in an interval of possible minimum point location and value by
 *  \a x and \a fx arrays. So, \a x location and \a fx value has an uncertainity. \n
 *  @f[ x[0] < real_x < x[1] @f]
 *  @f[ fx[0] < real_fx < fx[1] @f]
 *  We take the average of minimum and maximum values of uncertainity to bypass this uncertainity by AVE macro.
 *  This is not a huge problem because the uncertainity interval is so small. In most cases, minimum and maximum
 *  values of uncertainity are equal.\n\n
 *  <b>unexpected Cases</b>\n\n
 *  <b>Case 1: The three points are collinear and \f$ x[0] < x[1] < x[2] \f$ and \f$ fx[0] = fx[1] = fx[2] \f$ </b>\n\n
 *  If \ref Min_Search_Parabolic_Interpolation function returns -1 (Three points are collinear) but
 *  \f[ x[0] < x[1] < x[2] \f]
 *  \f[ fx[0] = fx[1] = fx[2] \f]
 *  Then, there is no error for us, because \a fx[1] value on \a x[1] point is the value that we are looking for.\n\n
 *
 *  <b>Case 2: The initial points failed to satisfy the condition that
 *             \f$ x1 < x2 < x3 \f$ and \f$ fx1 > fx2 < fx3 \f$ but
 *             \f$ x[0] < x[1] < x[2] \f$ and \f$ fx[0] = fx[1] = fx[2] \f$ </b> .\n\n
 *  If \ref Min_Search_Parabolic_Interpolation function returns -3, That is to say, if initial points failed to satisfy the
 *  the conditions below
 *  \f[ x1 < x2 < x3 \f]
 *  \f[ fx1 > fx2 < fx3 \f]
 *  but result values satisfies the conditions below:
 *  \f[ x[0] < x[1] < x[2] \f]
 *  \f[ fx[0] = fx[1] = fx[2] \f]
 * Then, The result values are valid for us. Becaue in some cases, we can come across some situations like below:
 *  \f[ x1 < x2 < x3 \f] and
 *  \f[ fx1 = fx2 < fx3 \f] or
 *  \f[ fx1 > fx2 = fx3 \f] or
 *  \f[ fx1 = fx2 = fx3 \f]
 *
 *  Then, there is no error for us, because \a fx[1] value on \a x[1] point is the value that we are looking for.\n\n
 *
 *  <b>Reduced Sample of \ref min_search function:</b>
 *  @code
 *  int min_search(double* fx, double* x, double* real_fx, double* real_x,const char direction) {
 *    int err=0;
 *    calc_Parabol_coef(fx[0], fx[1], fx[2]);
 *    err = Min_Search_Parabolic_Interpolation(f, &x[0], &fx[0], tolerance);
 *      if (err < 0 ) {
 *        if(err == -1) if( COND1(x) && ARE_EQUAL(fx) ) err=0;
 *        if(err == -3) if( COND1(x) && ARE_EQUAL(fx) ) err = 0;
 *        if(err<0) return err;
 *      }
 *    *real_x = AVE(x);
 *    *real_fx = AVE(fx);
 *    return err;
 *  }
 *  @endcode
 *  @param[in] *fx This is pointer of a a double array has 3 elements. ( \a fx[0], \a fx[1], \a fx[2] )\n
 *                 These elements holds the values of each point represented by <b> x location array </b> below.
 *  @param[in] *x  This is pointer of a double array has 3 elements. ( \a x[0], \a x[1], \a  x[2] )\n
 *                 These elements holds the coordinate values of \a fx values on the same axis(x or y)\n\n
 *                 <b>For Example</b>
 *                 \code
 *                  double x[3] = { -2.5, 0.0, 2.5 };
 *                  double fx[3] = { 67.0, 38.0, 76.0 };
 *                  //fx[0] value is on x[0] location
 *                  //fx[1] value is on x[1] location
 *                  //fx[2] value is on x[2] location
 *                 \endcode
 *  @param[out] *real_fx This is minimum value that found by the min_search function.
 *                       Function applies bi-parabolic interpolation method to the \a x and \a fx arrays and
 *                       finds a minimum \fx value between them. This value is assigned to \a real_fx variable.
 *                       This value should be between \a fx[0] and \a fx[2] and in mot cases must be smaller than \a fx[1] .
 *  @param[out] *real_x This is minimum location of \a real_fx variable. This value should be between \a x[0] and \a x[2].
 *  @param[out] direction This is 'x' or 'y' (Axis direction)\n
 *                        This variable is used ony for logging purposes and does not do nothing mch more than printing
 *                        operated axis name to the log files.
 *  @return Status of search procedure result.
 * @retval  0   Success
 * @retval -1   The three points are collinear.
 * @retval -2   The extremum of the parabola through the three points is a maximum.
 * @retval -3   The initial points failed to satisfy the condition that
 *              \f$ x1 < x2 < x3 \f$ and \f$ fx1 > fx2 < fx3 \f$ . */
int min_search(double* fx, double* x, double* real_fx, double* real_x,const char direction) {
	int err=0; /*Error handler*/
	log_deb("\nFor %c-axis:",direction);
	log_deb("%c={%.1f, %.1f, %.1f} and f%c= {%.0f, %.0f, %.0f}",direction,x[0],x[1],x[2],direction,fx[0], fx[1], fx[2]);
	//if( !ARE_EQUAL(fx) ) {
		f_parabola_calc_coefficients(fx[0], fx[1], fx[2]);
		err = Min_Search_Parabolic_Interpolation(f_parabola, &x[0], &fx[0], tolerance);
		if (err < 0 ) {
			if(err == -1){
				if( COND1(x) && ARE_EQUAL(fx) ) {
					err=0;
				}else{
					log_deb("\nError in input data %c-direction: err = %d", direction, err);
					log_deb("The three points are collinear.");
					if(fx[0] == fx[1]) log_deb("f%c[0] = f%c[1]",direction,direction);
					if(fx[1] == fx[2]) log_deb("f%c[1] = f%c[2]",direction,direction);
					if(fx[0] == fx[2]) log_deb("f%c[0] = f%c[2]",direction,direction);
					log_deb("f%c[0]=%.6f, f%c[1]=%.6f, f%c[2]=%.6f", direction, fx[0], direction, fx[1], direction, fx[2]);
				}
			}
			if(err == -2) {
				log_deb("\nError in input data %c-direction: err = %d", direction, err);
				log_deb("The extremum of the parabola through the three points is a maximum.");
			}
			if(err == -3) {
				if( COND1(x) && ARE_EQUAL(fx) ) {
						err = 0;
				}else{
					log_deb("\nError in input data %c-direction: err = %d", direction, err);
					log_deb(" The initial points failed to satisfy the condition that:");
					log_deb("%c1 < %c2 < %c3 and f%c1 > f%c2 < f%c3",direction, direction, direction, direction, direction, direction);
					log_deb("%12.8e < %12.8e < %12.8e and %12.8e > %12.8e < %12.8e",x[0],x[1],x[2],fx[0], fx[1], fx[2]);
				}
			}
			if(err<0) return err;
		}
	/*}*/
	/*We found the x-location of the real min. point*/
	*real_x = AVE(x);
	*real_fx = AVE(fx);
	log_deb("After parabolic-interpolation:");
	log_deb("%c={%12.8e, %12.8e, %12.8e} and f%c= {%.10f, %.10f, %.10f}", direction, x[0], x[1], x[2],direction, fx[0], fx[1], fx[2]);
	log_deb("%c-axis local minimum is %c=%12.8e f%c=%.10f",direction, direction, *real_x, direction, *real_fx);
	return err;
}

/*! @brief Search for the min point by parabolic interpolation and its value by biparabolic interpolation
 *
 *  If the hgt is a min point then use this function to obtain the real min point and location.
 *  This function applies the parabolic and biparabolic methods to find the results.
 *
 *  @param[in]  *h This is unpacked geoptantial height value on a specified grifd point in the NC file.
 *  @param[in]  scaleFactor Scale Factor for Geoptantial Height value conversion
 *  @param[in]  addOffset Add offset for Geoptantial Height value conversion
 *  @param[out] *x_loc is pointer to the x-location of the real min. value calculated by parabolic interpolation.
 *  @param[out] *y_loc is pointer to the y-location of the real min. value calculated by parabolic interpolation.
 *  @param[out] *val is the real min. value calculated by biparabolic interpolation method.
 *
 *  @return Status of search procedure
 *  @retval  0   Success
 *  @retval -1   The three points are collinear.
 *  @retval -2   The extremum of the parabola through the three points is a maximum.
 *  @retval -3   The initial points failed to satisfy the condition that
 *              \f$ x1 < x2 < x3 \f$ and \f$ fx1 > fx2 < fx3 \f$ . */
int hgt_find_real_min_location(short *h, float scaleFactor, float addOffset ,  double *x_loc, double *y_loc, double *val) {

	int err=0; /*Error handler*/
	double real_x=0.0;
	double real_y=0.0;/*0 means real min. point is exactly equal to middle value.*/
	double real_fx=0.0;
	//double real_y_val=0.0;
	double real_fxy = 0.0;
	/*
	 * We need to convert the height value's type from short to double and
	 * convert encrypted height values to real height values and store them in fx array.
	 * I created 3 array to store the 3 line with 3 elements.
	 */
	double Line[3][3];/*3x3 Node Matrix this will be modified by Min_Search_Parabolic_Interpolation routine.*/
	/* The zij matrix will be used to calculate the height values
	 * of two points in the positive and negative y-direction.*/
	double zij[3][3];
	double fx[3];
	double x[3] = { -2.5, 0.0, 2.5 };
	double y[3] = { -2.5, 0.0, 2.5 };
	register int i,j;
	for(j=-1;j<=1;j++) {
		for(i= -1;i <= 1;i++) {
			/*Convert encrypted height value to real height value*/
			Line[j+1][i+1] = UNPACK(VAL(h,i,j));
			zij[j+1][i+1] = Line[j+1][i+1];
		}
	}
	fx[0] = Line[1][0]; fx[1] = Line[1][1]; fx[2] = Line[1][2];

	log_deb("\n\tx0\tx1\tx2");
	log_deb("y0\t%.0f\t%.0f\t%.0f",Line[0][0], Line[0][1], Line[0][2]);
	log_deb("y1\t%.0f\t%.0f\t%.0f",Line[1][0], Line[1][1], Line[1][2]);
	log_deb("y2\t%.0f\t%.0f\t%.0f",Line[2][0], Line[2][1], Line[2][2]);

	err = min_search(fx, x, &real_fx, &real_x,'x');
	if(err<0) return err;

	/* We need to reset axis values to default values because
	 * Min_Search_Parabolic_Interpolation function changes them*/
	RESET_AXIS(x);
	double fy[3];
	/*
	 * Now, we need to calculate the height values of two points in the positive
	 * and negative y-direction which are y[0] = -2.5 and y[2]=2.5 points.
	*/
	fy[0] = Find_Zxy_Biparabolic_Interpolation(real_x,-2.5,&x[0],&y[0],&zij[0][0]);
	fy[2] = Find_Zxy_Biparabolic_Interpolation(real_x,2.5,&x[0],&y[0],&zij[0][0]);
	fy[1] = real_fx;

	err = min_search(fy, y, &real_fxy, &real_y,'y');
	if(err<0) return err;

	log_deb("local minimum is x=\"%.12f\" y=\"%.12f\"  fxy=\"%f\"", real_x, real_y, real_fxy);

	*x_loc = real_x;
	*y_loc = real_y;
	*val = real_fxy;

	log_deb("\nSUCCESS");

	return 0;
}

/*! @brief Search the minimum value
 *
 *  If the hgt is a min point then use this function to obtain the real min point and location.
 *  This function applies the parabolic and biparabolic methods to find the results.
 *
 *  @param[in]  *h is pointer to the min. height value will be used to find the real min. value and location
 *  @param[in]  scaleFactor Scale Factor for Geoptantial Height value conversion
 *  @param[in]  addOffset Add offset for Geoptantial Height value conversion
 *  @param[out] *x_loc is pointer to the x-location of the real min. value calculated by parabolic interpolation.
 *  @param[out] *y_loc is pointer to the y-location of the real min. value calculated by parabolic interpolation.
 *  @param[out] *val is the real min. value calculated by biparabolic method.
 *
 *  @return Status of search procedure
 *  @retval  0   Success
 *  @retval -1   The three points are collinear.
 *  @retval -2   The extremum of the parabola through the three points is a maximum.
 *  @retval -3   The initial points failed to satisfy the condition that
 *              \f$ x1 < x2 < x3 \f$ and \f$ fx1 > fx2 < fx3 \f$ . */
int hgt_find_real_min_location2(short *h, float scaleFactor, float addOffset ,  double *x_loc, double *y_loc, double *val) {

	int err=0; /*Error handler*/
	double real_x=0.0;
	double real_y=0.0;/*0 means real min. point is exactly equal to middle value.*/
	double real_fx=0.0;
	//double real_y_val=0.0;
	double real_fxy = 0.0;
	/*
	 * We need to convert the height value's type from short to double and
	 * convert encrypted height values to real height values and store them in fx array.
	 * I created 3 array to store the 3 line with 3 elements.
	 */
	double Line[3][3];/*3x3 Node Matrix this will be modified by Min_Search_Parabolic_Interpolation routine.*/
	/* The zij matrix will be used to calculate the height values
	 * of two points in the positive and negative y-direction.*/
	double zij[3][3];
	double fx[3];
	double x[3] = { -2.5, 0.0, 2.5 };
	double y[3] = { -2.5, 0.0, 2.5 };
	register int i,j;
	for(j=-1;j<=1;j++) {
		for(i= -1;i <= 1;i++) {
			/*Convert encrypted height value to real height value*/
			Line[j+1][i+1] = UNPACK(VAL(h,i,j));
			zij[j+1][i+1] = Line[j+1][i+1];
		}
	}
	fx[0] = Line[1][0]; fx[1] = Line[1][1]; fx[2] = Line[1][2];

	log_deb("\n\tx0\tx1\tx2");
	log_deb("y0\t%.0f\t%.0f\t%.0f",Line[0][0], Line[0][1], Line[0][2]);
	log_deb("y1\t%.0f\t%.0f\t%.0f",Line[1][0], Line[1][1], Line[1][2]);
	log_deb("y2\t%.0f\t%.0f\t%.0f",Line[2][0], Line[2][1], Line[2][2]);
	log_deb("\nFor x-axis:");
	log_deb("x={%.1f, %.1f, %.1f} and fx= {%.0f, %.0f, %.0f}",x[0],x[1],x[2],fx[0], fx[1], fx[2]);

	if( !ARE_EQUAL(fx) ) {
		/*
		 * Now, we will apply a parabolic interpolation on the x-direction includes
		 * fx[1] middle point in the heights array. Note that y is still 0.
		 */
		f_parabola_calc_coefficients(fx[0], fx[1], fx[2]);/*Firstly, calculate the coefficients of parabol*/
		err = Min_Search_Parabolic_Interpolation(f_parabola, &x[0], &fx[0], tolerance);
		if (err < 0 ) {
			if(err == -1){
				if( COND1(y) && ARE_EQUAL(fx) ) {
					err=0;
				}else{
					log_deb("\nError in input data x-direction: err = %d",err);
					log_deb("The three points are collinear.");
					if(fx[0] == fx[1]) log_deb("fx[0] = fx[1]");
					if(fx[1] == fx[2]) log_deb("fx[1] = fx[2]");
					if(fx[0] == fx[2]) log_deb("fx[0] = fx[2]");
					log_deb("fx[0]=%.6f, fx[1]=%.6f, fx[2]=%.6f",fx[0], fx[1], fx[2]);
				}
			}
			if(err == -2) {
				log_deb("\nError in input data x-direction: err = %d",err);
				log_deb("The extremum of the parabola through the three points is a maximum.");
			}
			if(err == -3) {
				if( COND1(x) && ARE_EQUAL(fx) ) {
						err = 0;
				}else{
					log_deb("\nError in input data x-direction: err = %d",err);
					log_deb(" The initial points failed to satisfy the condition that:");
					log_deb("x1 < x2 < x3 and fx1 > fx2 < fx3");
					log_deb("%12.8e < %12.8e < %12.8e and %12.8e > %12.8e < %12.8e",x[0],x[1],x[2],fx[0], fx[1], fx[2]);
				}
			}
			if(err<0) return err;
		}
	}
	//
	/*We found the x-location of the real min. point*/
	real_x = AVE(x);
	real_fx = AVE(fx);
	//double abs_real_x = ABS(real_x);
	//if(tolerance > abs_real_x) real_x = 0.0;
	//if(real_x==0) real_fx = zij[1][1];
	log_deb("After parabolic-interpolation:");
	log_deb("x={%12.8e, %12.8e, %12.8e} and fx= {%.10f, %.10f, %.10f}",x[0],x[1],x[2],fx[0], fx[1], fx[2]);
	log_deb("x-axis local minimum is x=%12.8e fx=%.10f", real_x,real_fx);
	/* We need to reset axis values to default values because
	 * Min_Search_Parabolic_Interpolation function changes them*/
	RESET_AXIS(x);
	double fy[3];
	/*
	 * Now, we need to calculate the height values of two points in the positive
	 * and negative y-direction which are y[0] = -2.5 and y[2]=2.5 points.
	*/
	fy[0] = Find_Zxy_Biparabolic_Interpolation(real_x,-2.5,&x[0],&y[0],&zij[0][0]);
	fy[2] = Find_Zxy_Biparabolic_Interpolation(real_x,2.5,&x[0],&y[0],&zij[0][0]);
	fy[1] = real_fx;

	log_deb("\nFor y-axis:");
	log_deb("y={%.1f, %.1f, %.1f} and fy= {%.10f, %.10f, %.10f}", y[0], y[1], y[2], fy[0] , fy[1] , fy[2] );
	if( !ARE_EQUAL(fy) ) {
		f_parabola_calc_coefficients(fy[0], fy[1], fy[2]);
		err = Min_Search_Parabolic_Interpolation(f_parabola, &y[0], &fy[0], tolerance);
		if (err < 0 ) {
		log_deb("\nError in input data y-direction: err = %d",err);
		if(err == -1) {
			if( COND1(y) && ARE_EQUAL(fy) ) {
				err=0;
			}else{
				log_deb("The three points are collinear.");
				if(fy[0] == fy[1]) log_deb("fy[0] = fy[1]");
				if(fy[1] == fy[2]) log_deb("fy[1] = fy[2]");
				if(fy[0] == fy[2]) log_deb("fy[0] = fy[2]");
				log_deb("fy[0]=%15.8e, fy[1]=%15.8e, fy[2]=%15.8e",fy[0], fy[1], fy[2]);
			}
		}
		if(err == -2) log_deb("The extremum of the parabola through the three points is a maximum.");
		if(err == -3) {
			if( COND1(y) && ARE_EQUAL(fy) ) {
					err = 0;
			}else{
				log_deb(" The initial points failed to satisfy the condition that:");
				log_deb("y1 < y2 < y3 and fy1 > fy2 < fy3");
				log_deb("%12.8e < %12.8e < %12.8e and %.10f > %.10f < %.10f",y[0],y[1],y[2],fy[0], fy[1], fy[2]);
			}
		}
		if(err<0) return err;
	}
	}
	//
	log_deb("After parabolic-interpolation:");
	log_deb("y={%12.8e, %12.8e, %12.8e} and fy= {%12.8e, %12.8e, %12.8e}", y[0], y[1], y[2], fy[0] , fy[1] , fy[2] );
	/*We found the x-location of the real min. point*/
	real_y = AVE(y);
	real_fxy = AVE(fy);
	log_deb("local minimum is x=\"%.12f\" y=\"%.12f\"  fxy=\"%f\"", real_x, real_y, real_fxy);
	*x_loc = real_x;
	*y_loc = real_y;
	*val = real_fxy;
	log_deb("\nSUCCESS");
	return 0;
}

//}//@}
